/**
 * Copyright (c) 2012, David Varbel
 * All rights reserved.
 * 
 * Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
 * 
 * 1) Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
 * 2) Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
package com.epistrata.neutraltheoryapplet;
/**
 * This code was inspired by Peter Abeles example at
 * http://code.google.com/p/efficient-java-matrix-library/wiki/PrincipleComponentAnalysisExample
 * 
 * Some sections of code were copied and slightly modified
 */

import org.ejml.alg.dense.decomposition.DecompositionFactory;
import org.ejml.alg.dense.decomposition.SingularValueDecomposition;
import org.ejml.data.DenseMatrix64F;
import org.ejml.ops.CommonOps;
import org.ejml.ops.SingularOps;

public class PrincipalComponentAnalysis {
	DenseMatrix64F original;
	DenseMatrix64F basis;
	double[] means;
	private int numComponents;
	
	public PrincipalComponentAnalysis(double[][] original, int numComponents) {
		if (numComponents > original[0].length) {
			throw new IllegalArgumentException("More components than dimensions in the input data.");
		}
		if (numComponents > original.length) {
			throw new IllegalArgumentException("Need at least as many samples as desired components in order to do PCA.");
		}
		this.original = new DenseMatrix64F(original);
		this.numComponents = numComponents;
		means = new double[original[0].length];
	}
	
	public void computeBasis() {
		//compute means for each dimension
		for (int i = 0; i < original.getNumRows(); i++) {
			for (int j = 0; j < original.getNumCols(); j++) {
				means[j] += original.get(i, j);
			}
		}
		for (int j = 0; j < means.length; j++) {
			means[j] /= original.getNumRows();
		}
		
		//normalize original data
		for (int i = 0; i < original.getNumRows(); i++) {
			for (int j = 0; j < original.getNumCols(); j++) {
				original.set(i, j, original.get(i, j) - means[j]);
			}
		}
		
		/*
		 * use Singular Value Decomposition instead of COV matrix per
		 * Peter Abeles example at http://code.google.com/p/efficient-java-matrix-library/wiki/PrincipleComponentAnalysisExample
		 */
		 SingularValueDecomposition<DenseMatrix64F> svd = 
             DecompositionFactory.svd(original.numRows, original.numCols, false, true, false);
	     if( !svd.decompose(original) ) {
	    	 throw new RuntimeException("SVD failed");
	     }
	     basis = svd.getV(true);
	     DenseMatrix64F W = svd.getW(null);

	     // Singular values are in an arbitrary order initially
	     SingularOps.descendingOrder(null, false, W, basis, true);
	
	     // strip off unneeded components and find the basis
	     basis.reshape(numComponents, original.numCols, true);
	}
	
	public double[] getBasisVector(int which) {
        if( which < 0 || which >= numComponents )
            throw new IllegalArgumentException("Invalid component");

        DenseMatrix64F v = new DenseMatrix64F(1,original.numCols);
        CommonOps.extract(basis,which,which+1,0,original.numCols,v,0,0);

        return v.data;
	}
	
	public double[] sampleToEigenSpace(double sample[]) {
		if (sample.length != original.numCols) {
			throw new IllegalArgumentException("Unexpected sample length");
		}
        DenseMatrix64F mean = DenseMatrix64F.wrap(original.getNumCols(), 1, this.means);

        DenseMatrix64F sampleMatrix = new DenseMatrix64F(original.getNumCols(), 1, true, sample);
        DenseMatrix64F result = new DenseMatrix64F(numComponents,1);

        //normalize
        CommonOps.sub(sampleMatrix, mean, sampleMatrix);

        //convert to eigenspace
        CommonOps.mult(basis, sampleMatrix, result);

        return result.data;
	}

}
